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ABSTRACT 


Magnetic torquing of spacecraft has been an important mechanism for attitude control since 
the earliest satellites were launched. Typically a magnetic control system has been used for 
precession/nutation damping for gravity-gradient stabilized satellites, momentum dumping for 
systems equipped with reaction wheels, or momentum-axis pointing for spinning and momentum- 
biased spacecraft. Although within the small satellite community there has always been interest in 
inexpensive, light-weight, and low-power attitude control systems, completely magnetic control 
systems have not been used for autonomous three-axis stabilized spacecraft due to the large 
computational requirements involved. As increasingly more powerful microprocessors have 
become available, this has become less of an impediment These facts have motivated 
consideration of the all-magnetic attitude control system presented here. 

The problem of controlling spacecraft attitude using only magnetic torquing is cast into the 
form of the Linear Quadratic Regulator (LQR), resulting in a linear feedback control law. Since the 
geomagnetic field along a satellite trajectory is not constant, the system equations are time varying. 
As a result the optimal feedback gains are time-varying. Orbit geometry is exploited to treat 
feedback gains as a function of position rather than time, making feasible the onboard solution of 
the optimal control problem. In simulations performed to date, the control laws have shown 
themselves to be fairly robust and a good candidate for an onboard attitude control system. 
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INTRODUCTION 

Magnetic torquing has been used for spacecraft attitude control since the launch of the 
earliest satellites. There are many operational spacecraft which use magnetics for 
precession/nutation damping, momentum dumping, and large angle maneuvers [1-4]. Also, there 
are existing three-axis control systems which use magnetic torquing to maintain spin axis 
orientation of a pitch reaction wheel [5]. Use of magnetics has recently been suggested for 
libration damping and arbitrary yaw angle control of a gravity gradient stabilized satellite [7]. 

Complete three-axis attitude control has not been used because of the large computational 
requirements involved. As increasingly more powerful microprocessors have become available, 
computation has become less of an impediment. This has motivated consideration of the all- 
magnetic attitude control system presented here. 

We consider a rigid spacecraft without a gravity gradient boom. For the present 
discussion, we assume a perfect knowledge of the spacecraft attitude is available. Finally, we 
assume the spacecraft is in a near circular orbit and has a nadir pointing nominal attitude. 

The geomagnetic field, while essentially constant in an earth-fixed reference frame, is time- 
varying in the nominal body reference frame. This introduces a time-varying element into the 
linearized system of equations. Because the magnetic torque vector is constrained to always lie 
perpendicular to the local geomagnetic field vector, the system appears uncontrollable if fixed at 
any instant in time. (When considered over time, it is completely controllable.) This does, 
however, limit the achievable closed loop transient reaction speed to die order of magnitude of the 
magnetic field time variation. For a nadir pointing satellite in a polar orbit, this is about two cycles 
per orbit. The reaction speed should decrease with decreasing orbit inclination because the 
geomagnetic field variation decreases in magnitude. 


FORMULATION OF THE EQUATIONS OF MOTION 

Nomenclature 

LQR Linear Quadratic Regulator 

P, R, Y Pitch, Roll, Yaw (2-1-3 Euler angles) 

<£>o Mean angular rate of the local vertical reference frame. 

Ri(Q) Rotation matrix for an angle 0 about the i-axis. 

Cb<— s Rotation matrix from inertial (space) to actual body coordinates. 

Cn<_ s Rotation matrix from inertial (space) to nominal body coordinates. 

Cbc_n Rotation matrix from nominal body to actual body coordinates. 

H 2 B Angular velocity of spacecraft with respect to inertial frame, represented 

in the body frame. 


ifi 


Total external torque acting on spacecraft body, in body frame. 
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Idist.B 


External disturbance torque, excluding gravity gradient torque, 
represented in the body frame. 

ic 0n> B Control torque acting on spacecraft body, represented in body frame. 

R s Magnitude of spacecraft position vector in Earth centered coordinates. 

A 

Rb Unit vector from center of Earth to spacecraft, represented in the body 

frame. 

J Spacecraft inertia tensor, represented in the body frame. 

Ji {J1.J2.J3} i* diagonal element of inertia matrix J. 

de(t) Control dipole moment of spacecraft, represented in the body frame. 

Hn Local geomagnetic field intensity, represented in the body frame. 

Hjsf(t) Local geomagnetic field intensity, represented in the nominal body 

frame. 

| 1 Geocentric gravitational constant, GMe = |i = 3.986 x 10 14 m3/s 2 . 


hi(t) {h 1 (t),h 2 (t),h 3 (t) } i* component of Hn( 0- 
di(t) { dj (t) ,d 2 (t) ,d 3 (t) } i* component of 4 b (0- 


S(0 

F 

G(t) 

J(dB(0) 

Sf 

A 

B 

S(t) 

K(t) 

K(p,X) 

a 

e 

i 


State vector [r, P, Y, R, P, y] 7 . 

System matrix of linearized equations. 

Input coupling matrix of linearized equations. 
Quadratic cost functional. 

Final time state cost matrix. 

State cost weighting matrix. 

Control cost weighting matrix. 

Solution of the matrix Riccati equation. 
Feedback gain matrix, parameterized by time. 
Feedback gain matrix, parameterized by (P,X). 
Orbit semi-major axis 
Orbit eccentricity. 

Orbit inclination. 



£2 Argument of the ascending node, 

to Argument of perigee. 

®E Angular rate of the Earth * 2k radians/24 hours. 

X Earth fixed longitude of the ascending node. 

P Argument of latitude (defined in Figure 4). 

FORMULATION OF THE EQUATIONS OF MOTION 

For the problem at hand, we have considered the rigid-body equations of motion of a nadir- 
pointing spacecraft in a circular orbit. Referring to Figure 1, where (S) subscripts indicate inertial 
frame, (N) indicates nominal body frame, and (B) indicates actual body frame, we can represent 
attitude using 2-1-3 Euler angles as 

Cb<— s = R 3 (Y) Ri(R) R 2 ((0 o t + P) (1) 

C B <_ s = R 3 (Y) Ri(R) R 2 (P) R 2 (co 0 t) (2) 

where Rj(0) is a rotation about axis i through angle 0. Since the desired trajectory has Y, R, and P 
all zero, we can define an intermediate "nominal" attitude (to be used later) as 


Cn<_s = R2(co Q t) , 


(3) 


which gives 


C B < — n = R 3 (Y) Ri(R) R 2 (P) . 


(4) 


The angular velocity of the body frame with respect to the inertial frame can be expressed as 


(Ob = Y 

i 

0 o 

1 

+ R 3 (Y) r 

' l ' 
0 

+ R 3 (Y)R!(R) (co 0 + P) 

o' 

1 


1 J 


. 0. 


- 0. 


To first order in angles and rates this is 


tito = 


R + co 0 Y 
C0 o + P 
L Y - co 0 R 


( 6 ) 


To obtain the first-order approximations to Euler's equations we write 

Ib = ~ (■ I £te) + tilB x J m , 
dt 

which, assuming a diagonal inertia matrix, gives 
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( 8 ) 


Ib ~ 


Ji R + ojo(Ji + J 3 - J 2 ) Y + tflg(J 2 - J 3 ) R 
J 2 P 

L J 3 Y - co 0 (Ji + J 3 - J 2 ) R + 0 )q(J 2 - Ji) Y 


The applied torque ib consists of the gravity gradient, control, and disturbance torques, and can be 
represented as 


IB — Idist.B Icon,B ^ [^B ^ (J £b)] > 
Rs 


(9) 


where R s is the geocentric radius, and Rb is the unit vector from the earth to the spacecraft, 
represented in the body frame. Since Rb corresponds to the z-axis in the nominal reference frame, 
we can use the first-order approximation 


Rb = 


-P 
R 
1 J 


( 10 ) 


giving 


3 H 


Ib — Idist,B + Icon,B " Q 

Rf 


(J 2 - J 3 ) R 
(Ji - J 3 ) P 
0 


(ID 


Finally, we examine the control torque icon.B- This torque is obtained from the 
commanded dipole jIbW as 

Icon.B = dfi(t) x Hb (12) 

where Hb is the ambient field intensity represented in the body frame. We also have 

HB = C B <_NH N (t) = H N (t) (13) 

where the explicit time dependence applied to Hn( 0 emphasizes that its value is attitude 
independent. We would like to replace Hb in equation (12) with Hn( 0. This is valid as a first- 
order approximation because we will later be assigning a control law for dB(t) which is linear in 
angles and rates. 

Combining equations (8), (11), and (12) with the approximation (13), we obtain the 
complete linearized equations of motion as follows: 


i(t) = • 



R 



P 


A 

Y 


dt 

R 



P 



_ Y - 



0 

0 

0 

F41 

0 

0 


0 

0 

0 

0 

F52 

0 


0 

0 

0 

0 

0 

F63 


1 

0 

0 

0 

0 

F64 


0 

1 

0 

0 

0 

0 


0 

0 

1 

F46 

0 

0 



R 


P 


Y 


R 


P 


_ Y - 
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0 

0 

0 

0 

G51(t) 

G61(t) 


0 

0 

0 

G42(t) 

0 

G62(t) 


0 

0 

0 

G43(t) 

G53(t) 

0 


di(t) 

d 2 (t) 

d 3 (t)J 


+ Idist.B 


where 


j i 1 R U 

G51(t) = -h 3 (t)J 2 1 

F52 = - — ( Jl ' j3 j 
R 3 \ h 1 

G61(t) = h 2 (t) J3 1 

F63 =^l 

G42(t) = h 3 (t) Jj 1 

F64-co 0 (Ji±M2.j 

G62(t) = -hj(t) J3 1 

F46 - -Wo[ Jl+ j 3 ' J2 ) 

G43(t) = -h 2 (t)JT 1 1 


G53(t) = hi(t) J2 1 


where 


hi(t), h2(t), h3(t) are the components of Hn (0 , 
di(t), d2(t), d3(t) are the components of 4 b (0 . 
or, in more compact notation, 

= Fi + G(t) dB(t) + Idist.B 


( 14 ) 


APPLICATION OF THE LINEAR QUADRATIC REGULATOR (LQR) 

Since it is desirable to maintain R, P, and Y as close to zero as possible over time, a 
reasonable cost function is the standard quadratic performance index [6], We wish to minimize 
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( 15 ) 


J(dB(t)) = lini x T (tf) Sf x(tf) +^- I [x T (t) A x(t) + d T e(t) B da(t)] dt 
tf ->~ 2 2 J i 0 

subject to the state equation 

*(t) = Fx(t) + G(t)d B (t) , 

where A is positive definite and B is non-negative-definite, assuming that idi st ,B = 0 in equation 
(14). This problem has a well-known solution given by 

d B (t) = -B- 1 G T (t)S(t)x(t) , (16) 

where the square matrix S(t) satisfies the Riccati differential equation 

S(t) = -S(t)F - F t S ( t) + S(t)G(t)B' , G T (t)S(t) - A (17) 

with the terminal condition S(tf) = Sf. This equation can be integrated backward in time to give the 
time-varying feedback gain matrix K(t) = -B~'GT(t)S(t). The dependence of S(t) on Sf becomes 
negligible as tf becomes large, a fact which can be exploited to give K(t) as if we were actually 
solving the problem with tf -» ». That is, tf is chosen far enough ahead that its effect is not seen at 
the present, and Sf, which can be any appropriately dimensioned non-negative-definite symmetric 
matrix, is set to zero. 

Kalman [8] has shown for a time-varying linear system, that if controllability is uniform 
over time, the closed loop infinite-time-horizon LQR is exponentially stable. 


RESULTS OF NUMERICAL SIMULATIONS 

The algorithm outlined above was simulated on an AT-class personal computer. The 
Riccati equation was solved numericallly using a fourth-order Runge-Kutta integrator, with time 
scaling to keep the system numerically well conditioned. The geomagnetic field was modelled 
using an 8th-order spherical harmonic expansion, as described in Wertz [9]. Nonlinear equations 
for rigid body dynamics and gravity gradient torques were used in the plant model. The feedback 
control law, while derived for the linearized plant equations, was applied to this nonlinear plant. 
The orbit used was near-circular and near-polar. 

Cost function matrices A and B were chosen empirically to produce desirable performance. 
Attention, however, was restricted to matrices of the form 


ae 


A = 


% 


ae 


a© 


a w 


a<o 


, and B - 


b 


b 


b J 


where a e , a^, b > 0. 
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In all cases tested, the closed loop system was stable. A simple check of controller 
sensitivity was made by using a slightly different model for the closed loop simulation than was 
assumed in computing the feedback gain. Moments of inertia, orbit eccentricity, and perigee 
position were varied, and in each case the system remained stable, though performance was 
degraded slightly, as expected. 

The satellite mass properties, orbit, and disturbance parameters used in evaluations are 
given in Table 1. The disturbance torques used are based on a 2 m 2 area and a 30 cm center of 
pressure to center of mass offset. This is conservative for a satellite with a maximum moment of 

inertia of 27 kg-m2. Figures 2A through 2E show time histories of Euler angles for several 
configurations. 


APPLICATION OF GEOMETRIC PROPERTIES 

Direct implementation of the time-varying linear-feedback control law described above 
requires a numerical solution of the matrix Riccati differential equation. Since the equation is 
solved backward from a future time to compute feedback gains at the present time, the integration 
must be restarted repetitively, each starting at a point later in time, as illustrated in Figure 3. 

Obviously, integration backward in time is not ideal for real-time implementation. 
Fortunately, however, we can utilize the quasi-periodicity due to orbit geometry to advantage, 
making real-time implementation feasible. ' ' 

We wish to parameterize the feedback gain K(t) in terms of position rather than time. First, 
we make certain assumptions about the orbit: 


(1) 

Orbit has a known semi-major axis, 

a 

(2) 

Orbit is circular, 

e = 0 

(3) 

Orbit has known inclination. 

i 

(4) 

Node is inertially fixed, 

dO/dt = 0 

(5) 

Perigee is inertially fixed. 

dco/dt = 0 

(6) 

Orbital angular rate is known, 

©o 

(7) 

Earth rotational rate is known, 

COE 


These assumptions will typically be valid over several days, which is adequate for our purposes. 
Referring to Figure 4, it becomes simple to describe an entire satellite trajectory, r(t), in Earth fixed 
coordinates, in terms of only p(to) and X(to) at any fixed time 


l(0 = R3(«E(t-to) - X(to)) Ri(-i) 


R 3 (-p(to) - COo(t-to)) 


■ 1 ■ 

0 

- 0 . 


( 18 ) 


The insight to be gained here is that if the trajectory can be completely parameterized in terms of the 
ordered pair (P(to),X(t 0 )) then the geomagnetic field history, Hj\j(t) can be similarly expressed in 
terms of (p,X). Hence the feedback gain matrix K(t 0 ), which depends on H_n(0 from t 0 <t<°°, can 
be written as a function of (p(to),X(to)) only. This function, K(p,X), is not time varying. Thus, 
K(P,X) can be computed once and will remain valid permanently. 


A real orbit will not satisfy the assumptions exactly. We can, however, drop the 
assumption of a circular orbit; the expression in equation (18) then becomes much more tedious, 
but the conclusion is the same. That is, the trajectory can be completely characterized by (p,X). 
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When assumptions (4) and (5) are violated by a small amount, the effect is to add a slow 
time variation to K(P,X). In an implementation of this system, K(p,k) must be recomputed 
periodically. It is expected that recomputing once every three days will be adequate for a typical 
perigee drift rate of three degrees/day. This requires much less CPU time than the continuous gain 
computation of Figure 3. 

The feedback gain function K(p,X) has been computed for several cases. It appears to be 
reasonably well behaved and can be represented with a modest investment in computer memory. 
Interestingly, a secondary advantage of using (p,k) to parameterize K(*), rather than another 
position representation scheme, is that K(p,k) is periodic in both p and X. with a period of 2k. It 
may be possible to reduce data storage requirements by representing K(p,X) in terms of Fourier 
coefficients. 


FLIGHT SOFTWARE CONFIGURATION 

A suggested flight software configuration is given in Figure 5. The infrequent updating of 
the feedback gains is run as a background task, and the Real Time Control operation is the high 
priority task. 

Since both tasks require a model of the orbit, periodic orbit updates from a ground station 
will be necessary. It is believed that updates on the order of once every seven days should be 
adequate for most applications. 


CONCLUSIONS 

The feasibility of spacecraft attitude control using only magnetic torquing has been 
demonstrated. Although the closed-loop transient reaction speed possible with such a system is 
fundamentally limited, the attitude requirements of many missions appear attainable. An algorithm 
for flight computer implementation has been simulated, demonstrating the feasibility of using this 
system with a typical onboard microprocessor. The mechanical simplicity inherent in using 
electromagnetics only for control promises to make such a system both cost effective and 
mechanically reliable. 
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Table 1 


System Parameters Used in Test Cases of Table 2 


Moments of Inertia 



It 

27 

kg-m 2 


h 

27 

kg-m 2 


13 

10 

kg-m 2 

Orbital Rate 


CD 0 

1.0473E-3 

rad/s 

Orbit 

Semi-major axis 

a 

1.15537 

earth radii 

Eccentricity 

e 

0.004119 

degrees 

Inclination 

i 

89.547 

Right ascension of 


111.167 

degrees 

the ascending node 

Arg. of perigee 

COp 

71.0035 

degrees 

Node drift rate 

clQ^/dt 

-0.047522 

deg/day 

Per. drift rate 
Epoch day 

dcop/dt 

-3.004395 

130 

deg/day 

Time of day 
Year 


78194.5 

87 

sec 

Orbservation Time Span 

Simulations start at 1987, day 131, time = 100 sec. 



Disturbances Used in Performance Analysis 


Torque (constant in body frame) 

±55 dyne-cm in Roll 

±100 dyne-cm in Pitch 

±10 dyne-cm in Yaw 


Residual Magnetic Dipole (constant in body frame) 
200 pole-cm in body x-axis 

200 pole-cm in body y-axis 

200 pole-cm in body z-axis 




X S = Y s xZs 
Ys = orbit normal 

Zs = center of earth to S/C position at time t=0 


X N = Y n x Zn 

Yn = orbit normal P = Pitch 

Zn = local vertical (up) Y = Yaw 


Figure 1: Coordinate System and Small 2-1-3 Euler Angle Definitions 



Figure 2a. history of Euler angles for a case without any disturbance torques or a residual 
£lydlc£ ^ CS d ° n0t dCCay precisely to zero ^ause the orbit is not 
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Figure 2b: Case where disturbance torque is a constant +10 dyne-cm in roll, +100 dyne-cm in 
pitch, and +5 dyne-cm in yaw. A residual dipole of 200 i + 200 j + 200 k pole-cm is 
also assumed. 


[EULER ANGLES: DlBt. -55 1QO “IQ dyne-cml 




Figure 2c: Disturbance torque for this case is a constant -55 dyne-cm in roll, +100 dyne-cm in 

pitch, and -10 dyne-cm in yaw. A residual dipole of 200 i + 200 j + 200 k pole-cm is 
also assumed. 
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Figure 2d: distmtance torque assumed for this case is +55 dyne-cm in roll, + 1 00 dyne-cm in 

assumed dyne ' Cm ln yaw ‘ A residual dipole of 200 i + 200 j + 200 k poie-cm is 



o South Pole .« North Pole i South Pole « « North Pole - 


TIME (orb i ta) 


F,gurc 2e: is a constant -55 dyne-cm in roll, +100 dyne-cm in pitch, and -10 

TTw rlncMt i^ aW ^ resl< iaal dipole of 200 i + 200 j + 200 k pole-cm is also assumed, 
l ne closed loop system will not capture reliably at initial angles larger than 40°. 
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Figure 3: Iteration necessary when using backward Riccati equation integration in real time. 



Notes: 

(1) ^ = -coe, cog = 2k radians/24 hours 
dt 

(2) For near circular orbit, ^-= oo 0 , co 0 = orbital angular rate (rad/sec) 

dt 


Figure 4: Definition of position parameters |J and X. 
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Figure 5: Control Software Structure 









